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Abstract 

We simulate Quantum Chromodynamics (QCD) in four Euclidean 
dimensions with two (degenerate mass) flavors of dynamical quarks. 
The Dirac operator we use is the so-called chirally improved Dirac 
operator. We discuss the algorithm used for the simulation as well 
as the checks and some results on lattices up to size 8 4 for fermion 
mass parameters down to 0.1. This is the first attempt to introduce 
dynamical quarks with the chirally improved Dirac operator. 
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1 Introduction 



Spontaneous breaking of the chiral symmetry is one of the key issues in QCD. 
For a lattice study of this phenomenon it is desirable to have a formalism 
which maintains chiral symmetry as much as possible. Massless lattice Dirac 
operators obeying the Ginsparg- Wilson condition (GWC) [I], 

D j + D = 2D^RD, (1) 

(with a local operator R), provide the weakest form of violation of chiral 
symmetry: it is violated locally and restored in the continuum limit. Also, 
since the smallest quark masses are a few MeV we have to approach the chiral 
limit in realistic lattice simulations eventually. 

Fermion zero modes are related to topological properties of QCD and 
may be important for spontaneous chiral symmetry breaking in the chiral 
limit. It is therefore of particular importance for simulations with dynamical 
fermions to be able to approach this limit. At the moment only fermions 
defined through Dirac operators obeying the GWC (GW- fermions) seem to 
be suited for an approach to the chiral limit. An important property of such 
operators (|TJ) for R = 1/2 is that their eigenvalue spectrum lies on a unit 
circle centered at 1 in the complex plane. 

It is extremely costly to include the full fermion dynamics with such 
fermions as compared to the simpler Wilson or staggered fermions formula- 
tion, with presently available algorithms. In fact, there have been very few 
attempts to use GW-fermions in dynamical simulations [21 El IH El E] and all 
of them have been exploratory in nature. 

These studies have been concerned with overlap fermions |7j which are 
the only known exact GW-fermions. This property, however, also gives rise 
to additional problems due to the discontinuous development of the operator 
spectrum even for a continuously changing gauge field. With the overlap 
operator, zero modes appear or disappear instantaneously accompanying an 
overall change of the Dirac operator spectrum. 

Computationally more economic solutions are domain- wall fermions [S], 
fixed-point fermionsjHJ or chirally improved fermions [TO] . Even though these 
actions either need an extra dimension or they have considerably more terms, 
they are typically a factor of O(10) more expensive than simpler actions but 
still O(10) less expensive than overlap fermions. 
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1.1 Chirally improved fermions 

Here we will work with the chirally improved fermions. Quenched simulations 
for this action have demonstrated good chiral properties allowing for pion 
masses down to ~ 250 MeV. The ground state hadron spectrum has been 
determined, e.g., in and for excited hadrons in |12j . 

The chirally improved massless Dirac operator may be written as a trun- 
cated series of terms 

Sci = $(x)D C i(x,x + P)iP(x + P) , 

x P 
16 

D cl {x,x + P) = ^2r a c%U(x,x + P) . (2) 

a=l 

The sum in the action runs over path shapes P connecting x with x + P 
while the sum over a in Da runs over all elements of the Clifford algebra. 
U denotes the ordered product of link variables along this path. 
The massive operator we define as 

Dci(m) =m + Dd, (3) 

where m denotes the dimensionless quark mass , i.e., it is the valence quark 
mass m va i in the quenched simulation and agrees with the sea-quark mass 
m sea in the dynamical case. 

Plugging (J2J) into the GWC and truncating the system (number of co- 
efficients and equations) one obtains a set of algebraic relations for the co- 
efficients dp. The lattice symmetries, invariance under charge conjugation 
and parity as well as 75-hermiticity are respected but the series is truncated 
at path length 4 and only a subset of 19 coefficients has been considered. 
These coefficients depend implicitly on the gauge coupling and have to be 
re-determined at different values. The leading gauge coupling dependence is 
- similar to tadpole improvement - coded in two parameters z s and z v which 
multiply the gauge links in the formal expansion. Usually it is sufficient to 
take z s = z v and one adjusts this value such that the spectrum of Da, which 
approximately follows the GW-circle, passes through zero. Coefficients for 
the quenched simulation [TT] can be found in [TB"] . 

1 Actually, due to a trivial renormalization the correct mass value is m/(l + to/2) but 
for simplicity of notation we always refer to to here. 
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This is an important technical point: The coefficients of Dqi depend 
on both, the gauge coupling /3 and the sea-quark mass m sea and have to be 
determined by adjusting z s and solving the algebraic equations resulting from 
the GWC as discussed in [TU] . 

Our Dirac operator is always defined on one-step hypercubic (HYP) 
smeared [H] gauge configurations in order to reduce ultraviolet (UV) fluctu- 
ations. In that sense the definition of our Dqi includes the smearing step. 



1.2 Liischer-Weisz action with tadpole improvement 

Previous experience in quenched calculations showed that using the Liischer- 
Weisz action with coefficients from tadpole improved perturbation theory 
leads to nicer chiral properties for this Dirac operator JU|. I n particular the 
spectrum of the Dirac operator at small eigenvalues deviates less from the 
circular shape. For that reason we also use that gauge action here; it reads 

Slw = -A 22 3 Re tr ^p 1 ^ ~ ^ 2 3 Re tr Urc 

plaq re 

-&V^Re trC/ tb) (4) 
tb 6 

where U p \ m is the usual plaquette term, U Te are Wilson loops of rectangular 
2x1 shape and U t b denote loops of length 6 along edges of 3-cubes ( "twisted 
bent" or "twisted chair"). The coefficient (5\ is the independent gauge cou- 
pling and the other two coefficients and (3% are determined from tadpole- 
improved perturbation theory. They have to be calculated self-consistently 
[To] from 

U ° = Q Re tr ^P la ^) ' a = " 3 06839 lo S W) > ( 5 ) 

through 

p 2 = (i + 0.4805 a) , fa = ^0.03325 a . (6) 

Again, this determination should be done for each pair of couplings m sea ). 

We discuss here results for lattices up to size 8 4 and sea-quark masses 
down to 0.1; our emphasis lies on the method, although we do discuss effects 
of dynamical fermions on lattice spacing and propagators. Like other studies 
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for GW-fermions, our study also has an exploratory character hopefully on 
the way towards implementing the approach for large scale simulation. 



2 The updating algorithm 

We simulate QCD with 2 flavors of quarks with degenerate sea-quark mass 
m sea using the chirally improved Dirac operator. The action thus has the 
form 



where is the usual pseudo-fermion field [TH] . 

As mentioned, our Dirac operator includes HYP-smearing of the gauge 
configuration. This smearing procedure involves the projection of a general 
complex matrix into SU(3). This operation is not different iable and the use 
of the exact Hybrid Monte Carlo (HMC) method ^7] is therefore ruled out. 

The algorithm we implement can be thought of as a variation of standard 
HMC. To ensure detailed balance in HMC one introduces auxiliary momenta 
p (conjugate to U) and defines a Hamiltonian Ti by 



where S[<p, U] is the original action of the theory. The molecular dynamics is 
driven by the Hamiltonian equations of motion and the pseudo-fermion field 
<fr is held fixed throughout the molecular dynamics trajectory. The final step 
is an accept/reject step with the acceptance probability P acc given by 



The equilibrium distribution is determined entirely by the action used in 
the accept /reject step as long as the molecular dynamics trajectories are 
reversible. The molecular dynamics evolution does not necessarily have to 
be generated by the same action. Exploiting this freedom we use a two-step 
algorithm in which the first step consists of making a proposal according 
to some simple (computationally cheap) action development and the second 
step is the accept/reject step with the original action. 

For the molecular dynamics step we define our Hamiltonian by 



S[<j), U] = S LW + ^(Dcr Vsea^PcrVsea))^ , 



(7) 



n= ~tr P 2 + s[4>,u], 



(8) 






(10) 
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where Sample is a simpler, numerically cheaper action. Our acceptance prob- 
ability nevertheless is given by 



P, = mi n(l,gg^±##4) (11) 
l exp(-{p 2 + S[<#>, U]} M ) ) 

where S now is the full, original action. Thus we have an exact algorithm and 
need not worry about systematic biases. The central problem now is how to 
generate proposals for configurations efficiently, such that the path through 
configuration space is as close as possible to that of the original action, i.e., 
how to choose an efficient Sample . 

On our path to the finally chosen algorithm we tested various alternatives. 
An obvious first choice for Sample is the Luscher-Weisz gauge action and the 
Wilson Dirac operator. The parameters of the Dirac operator were chosen to 
correspond approximately to the plaquette value and sea-quark mass repre- 
sented by our values chosen for Dqi in the accept/reject step. We then also 
tried to replace the Wilson Dirac operator by a truncated Dqi including only 
terms up to length 2. 

Finally we turned off the fermionic part of the molecular dynamics equa- 
tions. At that point our Sample consisted only of the Luscher-Weisz gauge 
action. This turned out to be superior: it is faster and the inclusion of the 
fermionic parts did not significantly improve the final acceptance. 

The accept/reject step involves the ratio of determinants of the Dirac 
operator. This is approximated by the stochastic estimator method inherent 
in the pseudo-fermion formulation. It was pointed out ^Hl EE 120] that the 
noise in this stochastic estimation may introduce artificial barriers on the 
way through configuration space. Various methods to reduce the fluctuation 
of this estimate have been discussed. Following these ideas we introduced 
an additional UV- filter |2T| . The basic concept is to reduce the spread of 
the eigenvalues of the operator to be stochastically estimated. One defines a 
reduced matrix 

D r = Dexp(/(D)), (12) 

where f{D) is chosen to be a polynomial in D with coefficients such that the 
eigenvalues of D r are concentrated around z — 1 in the complex plane. In 
our acceptance step we need to compute 

det(£>t D) ncw , /n , det(D r t D r ) new 

- exp (-2 tr f(D) aew + 2 tr /(-D)oid) - — — f : — • (13) 



det( J DT J D) old det{DlD T ] 



old 
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Figure 1: We compare the eigenvalue spectrum of D C i (+) with that of the 
reduced operator Z?ci,r ( x ) f° r a typical configuration (4 4 lattice, fti = 7.4, 
m sea = 0.1); the latter are much closer to the value 1. 



A complication for extended Dirac operators is to compute the trace over 
polynomials of D. However for Dq\ it is relatively straightforward to compute 
at least tr(D) and tr(D^D) . We tried a polynomial 

f(D) = 7 tr(D) + 5tv(D^D) , (14) 

but found that introducing 5 did not affect our results significantly. We 
therefore made the simplest choice f(D) — 7 D with 7 = —0.477. In Fig. ^ 
we compare the eigenvalue spectrum of Dq\ with that of the reduced operator 
and find, indeed, that the reduced operator is closer to unity. 

Our finally chosen updating method thus used this UV-nlter combined 
with HYP-smearing and molecular dynamics equations using only the LW 
gauge action. Our trajectory lengths are 0.07 units of molecular dynamics 
time. In the actual implementation our trajectory consisted of a half step, a 
full step and a final half step in terms of the gauge field. 
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The acceptance itself is done in two steps. In the first step the value 
exp (— 27tr (D new — D \ d )) is calculated exactly along with the change in 
the kinetic energy and bosonic part of the action. If this is accepted, the 
second - more expensive - step is the stochastic estimation of the ratio of 
determinants of the reduced matrices. The acceptance in the first step is 
~ 23% and that of the second step is ~ 20%. Although this gives an overall 
acceptance rate of less than 5%, one has to keep in mind that neither the first 
step nor the proposal involve any inversion of the fermion matrix and both 
are therefore quite fast. The only time consuming step is the second. Thus 
the net efficiency is to be compared with an HMC close to ~ 20% acceptance 
rate, albeit with smaller trajectory length. 

Since the inversion of the chirally improved operator is the most expensive 
step in our algorithm we tried to increase its efficiency as much as possible. It 
is well-known that during the trajectory development in standard HMC 
the inversion of the Dirac operator can be speeded up by using the solution 
vector of the previous step as the initial guess for the current step. Our case is 
slightly different. We do not invert the Dirac operator during our trajectory 
development, but our trajectory lengths are smaller. Therefore we expect 
that the new gauge field configuration U new is not too far away from the 
starting field configuration U Q \d- Also we note that the pseudo-fermion field 
is generated by = -Dci(koid) where £ is a complex gaussian random 
vector. Thus we have £ = -Dci(^oid) -1 4> an d we use this vector £ as an initial 
guess for the inversion. Indeed this choice reduces the necessary number of 
matrix-vector multiplications by about 20%. 

Another speeding-up technique we use is to utilize the fact that we do not 
need to estimate the determinant ratio exactly. If r] is the random number 
with which we want to compare the ratio, then we need to check only if 
(—log//) is larger than the change in action or not [T2|. For the overlap 
operator this can be implemented as follows. The Metropolis accept/reject 
step compares the norm of solution vector x = Dov" 1 <P with (— log 77). At 
the n-th step of the bi-conjugate gradient routine, let the solution vector be 
x n and the residual vector r n = — D Qv x n . For the overlap operator one 
knows that 

< |£>ov _ Vl < -101 • (15) 



2 + m m 
Then it is straightforward to show that 

2 1 22 2 2 1 2 

m [2 + m) z m m z 
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Since Dqi satisfies the GW- relation (with R = 1/2) to a good approxima- 
tion, we assumed (and checked numerically) that its spectrum satisfies these 
bounds too. So at every step we only need to compute the upper and lower 
bounds on \x\ 2 to see if (— log?y) is inside that range or not. This reduces 
the number of matrix vector multiplications typically by a factor 3. For a 
test of our assumptions we check the accuracy of this method by computing 
the solution vector to an accuracy of 10 -12 randomly once in 20 updates on 
the average. These tests never failed in our study. 

We checked our numerics in several ways. Internal consistency of the 
program was checked by verifying the reversibility of the molecular dynam- 
ics trajectories. After a forward and backward trajectory the final energy was 
equal to the starting energy up to the precision of our calculations (double 
precision). As another check of the possible problem of using a noisy estima- 
tor we compared the stochastic estimator for the ratio of determinant with 
an exact evaluation (on the 4 4 lattices). We found that the resulting plaque- 
tte expectation value was equal within errors (less than 0.15%) whether we 
calculated the determinant exactly or estimated it using the pseudo-fermions. 

Further checks, also on 4 4 lattices, were to reproduce the quenched pla- 
quette values by turning off the fermions, reproduce dynamical Wilson results 
by replacing the chirally improved operator by the Wilson Dirac operator in 
the accept/reject step and reducing the chirally improved operator to the 
Wilson Dirac operator by changing the coefficients. 

We have discussed in the introduction that the Dqi parameters depend 
on the normalization parameter z s which is adjusted such that the massless 
operator has eigenvalues running through zero. This parameter is a function 
of fa and m. Also the LW-action parameters fa and fa are functions of fa and 
m. All of them, z s (fa, m), fa(fa, m) and fa(fa, m) have to be determined 
self-consistently by iterating the defining equations for Dqi and the LW- 
action. 

To determine fa and fa self-consistently due to (J3J we used a "moving 
average" of the plaquette, i.e., the average of the plaquette over a reasonably 
large interval of successive updates, and set it to Uq. The interval is shifted 
with new updates by dropping the oldest point and adding a new one. The 
moving average has less fluctuations compared to the original plaquette and 
in equilibrium it is practically identical to the plaquette average. Once equi- 
librium is reached we do not change the gauge couplings any more. The final 
numbers are given in Tabled 
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3 Results 



In the quenched case for the LW-action the lattice spacing values for various 
values of (3\ have been determined in [221 I2H] • For (3\ = 7.6 we have a lattice 
spacing a « 0.19(1) fm, corresponding to a lattice size of 1.5 fm for the 8 4 
lattices in the quenched case. 

Our final results are mainly from a simulation of the chirally improved 
Dirac operator on a 8 4 lattice with the tadpole improved Luscher-Weisz gauge 
action at fl\ = 7.6 and quark mass parameter m sea = 0.1. The measured 
plaquette value for this run (cf. Tabled} compares very well with the assumed 
plaquette value used for the determination of the LW-action (jlj). Unless 
explicitly stated otherwise, all results refer to this run; it is a sequence of 
120000 updates. With our average acceptance rate of ~ 5% per update 
(i.e. per accept/reject step) and individual trajectory lengths of 0.07 this 
corresponds to a total effective molecular dynamics time of 420, counting only 
accepted steps. We allowed 40000 updates for equilibration and then saved 
a configuration every 2000 updates. All our propagators and masses have 
been computed on 40 such configurations. In computing the propagators an 
additional complication, compared to the quenched case, is that the masses 
of the sea-quarks and the valence- quarks should agree; one therefore cannot 
use multi-mass solvers. 

The total run time for equilibration and configuration generation was 
about 2 weeks on a Linux cluster using 32 2.4 GHz Xeon CPU's. 

3.1 Equilibration 

Equilibration for the chirally improved operator is a rather slow process. 
Also on a 4 4 lattice we noticed that the number of matrix vector multiplica- 
tions required were very high during the initial part of the cold start. This 
made cold starts quite impracticable for larger lattices. On the 8 4 lattice we 
chose for our starting configurations quenched configurations with plaquette 
values significantly higher and lower than the expected dynamical equilib- 
rium value and let the two sequences converge. The plaquette history for 
such a process is shown in Fig. where we denote the plaquette value by 
P= | Re tr(C/ plaq ). 

Another important question is that of autocorrelation. Experience with 
exact HMC shows that very small step sizes lead to rather large autocorre- 
lation times. We believe that with our trajectory length of At = 0.07, the 
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Figure 2: Equilibration of the plaquette for lattice size 8 4 , Pi = 7.6 and quark 
mass m sea = 0.1. 



autocorrelation times are moderate. This is based on a quenched study of 
the autocorrelation time of the plaquette and we show the results in Fig. El 

From this figure we see that while the run with the trajectory length of 
1.024 falls fastest, the run with trajectory length 0.064 is not too different 
from it after 20 time units whereas the run with length 0.016 is still quite far 
away. Assuming that such a picture also holds for the dynamical case, we 
conclude that our autocorrelation length is only moderately larger than for 
standard HMC where one typically uses a trajectory length ~ 1. 

To get an idea of the autocorrelation time in our runs, we did a binned 
error analysis for the plaquette. We plot the result in Fig. 0] As can be 
seen clearly from the figure, the maximal error is obtained around a bin size 
of 2000, corresponding to an effective molecular dynamics time interval of 7. 
We take this to be an estimate of our autocorrelation time. 

On a 4 4 lattice we also studied the dependence of the plaquette expec- 
tation value on the quark mass. These studies were carried out for the 
Luscher-Weisz action at f3i = 7.4. The results are given in Table ^ together 
with the values for the simulation on the 8 4 -lattice. Since switching on dy- 
namical fermions should in leading order be equivalent to going to larger /3i 
in a quenched simulation, one expects that the plaquette value increases for 
decreasing sea-quark mass. This is indeed what we observe. 
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Figure 3: Equilibration of the plaquette from cold start as discussed in the 
text. The three different curves correspond to trajectory lengths of 1.024, 
0.064 and 0.016 in molecular dynamics time. 



3.2 Spectrum of the Dirac operator 

The chirally improved Dirac operator is an approximate solution of the 
Ginsparg- Wilson relation. In particular it has the property that the low 
lying spectrum is close to the Ginsparg- Wilson circle. In order to verify that 
this property was preserved also with dynamical quarks, we plot the first 
100 eigenvalues for three equilibrium configurations in Fig. comparing 
the quenched with the dynamical situation. Compared to the quenched case 
the spectrum in the dynamical case is closer to the circle, indicating smaller 
effective lattice spacing. 

Zero modes of the Dirac operator are related to instantons. For the 
chirally improved operator, in the quenched as well as the dynamical case 
with large quark masses, zero modes are observed. However for the mass 0.1 
of our 8 4 simulation we did not find any zero modes after equilibration. 

In order to understand this feature, which seems to contradict other find- 
ings El we also ran simulations with the Wilson Dirac operator on a 4 4 
lattice. Here we did see would-be zero modes (i.e. reasonably small real 
eigenvalues) quite frequently. This may be due to the fact that the fluctua- 
tions of the Wilson Dirac operator seem to be much larger than the chirally 
improved operator. In fact with our algorithm we were not able to simulate 
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Figure 4: Variation of the error for the plaquette with bin size allowing for 
an estimate of the autocorrelation time (8 4 lattice, f3\ = 7.6 for the Liischer- 
Weisz gauge action.) 



Table 1: Measured plaquette values for different masses; Uq denotes the as- 
sumed plaquette value. 



V 


Pi 


^sea 


0.1 


0.5 


2.0 


10.0 


oo 


4 4 


7.4 


^plaq 


0.608(1) 


0.604(1) 


0.591(1) 


0.579(1) 


0.556(2) 






«0 


0.606 


0.601 


0.591 


0.582 


0.556 


8 4 


7.6 


^plaq 


0.6202(2) 








0.5824(1) 






U 


0.62 








0.5825 



the Wilson Dirac operator on the larger 8 4 lattice. The fluctuations of the 
fermionic determinant were far too large for any reasonable trajectory length. 

As a further check we also replaced the stochastic estimator for the deter- 
minant by the exact evaluation for the small lattice runs. The tunneling be- 
havior did not change. Even when we started with a topologically non-trivial 
quenched configuration (i.e., zero modes of Dqi for the quenched gauge con- 
figuration) the configurations quickly tunneled to the zero-topological charge 
sector while performing the HMC updates. We conclude that there is no 
obvious tunneling barrier in our algorithm but that the trivial zero mode 
sector is natural for our choice of couplings. 
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Figure 5: Spectra of the chirally improved operator on a 8 4 lattice. (3\ for 
the Liischer-Weisz gauge action was 7.6 and the quark mass m = 0.1. These 
spectra are for three typical, randomly selected quenched (diamonds) and 
dynamical (plus) configurations. 



Recent results on the dynamical overlap fermions do report seeing zero 
modes [SJ E] . However, the plaquette values quoted in [5] are quite different 
from ours; indications are that those runs are effectively at smaller (3 than 
ours. That coupling leads to a larger physical volume (and lower tempera- 
ture) and it is not obvious that the results can be compared. The results of 
[S| extend to similar parameter values as ours and tunneling to sectors with 
one zero mode were observed. We work at a slightly larger gauge coupling, 
and, as argued below, we may be deeper in the deconfined phase. More 
statistics and fine tuning of the parameters will be necessary to resolve this 
discrepancy. 

3.3 Propagators 

One of the primary goals of lattice simulations is to reproduce the known 
hadron mass spectrum. Although the lattice size is definitely too small to 
identify asymptotic states, we still have computed pion and rho meson prop- 
agators with our dynamical chirally improved configurations for a crude, 
preliminary check. The propagators were computed using point sources and 
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Figure 6: Normalized pion (full circles) and rho (open circles) correlators 
on a 8 4 lattice. We compare the quenched results (left-hand plot) with the 
dynamical ones (at m sea = 0.1, right-hand plot), both at Pi = 7.6 and 
valence-quark mass m val = 0.1. The curves represent the cosh-fits to the 
three central points. 



sinks. Our expressions for the pion and the rho correlators are given by 

CV(0,t) = J]tr(7 5J D c r 1 (x,t:0,0)75^cr 1 (0,0:x,t)) , (17) 

C p (0,t) = tx( li D c r 1 {x,t:0,0) li D cl - l {0,0:x,t)) . (18) 

£,1=1,2,3 

The correlation functions are shown in Fig. El for the dynamical (m sea = 
rrivai = 0.1) as well as for the quenched (m va i = 0.1) case. The error bars are 
naively determined without auto-correlation analysis. 

Some conclusions can be drawn comparing the propagators of the dy- 
namical with those of the quenched case. Table 121 shows the results for the 
"meson masses" m of our cosh(m (t— 7ix/2))-fits in the symmetry region. Let 
us denote the dimensionless masses by m q (for quenched) and m d (for dy- 
namical) and use the values from the smaller fit range (with better x 2 1 d.o.f). 
We find m^/ml m 1.52 and m d p /m q p ~ 1.04. There can be two reasons for 
this. Either the lattice spacing has increased in the dynamical case or it 
has decreased so much that we are at much smaller physical time extent 
and therefore cannot observe the asymptotic decay of the correlators. Our 
spectrum, as discussed in the previous section, suggests that the second ex- 
planation is more plausible. Thus we cannot use these values to derive the 
lattice spacing. 
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Table 2: Fitted mass values in lattice units for lattice size 8 4 , (3\ = 7.6 and 
valence-quark mass m va i = 0.1. 







Pion 


Rho 


simulation 


range of A t 


"V X 2 /d-o.f 


m p x 2 /d.o.f 


dynamical 


2- 6 

3- 5 


1.18 (3) 10.54 
0.97 (1) 0.02 


1.33 (4) 10.62 
1.06 (2) 0.14 


quenched 


2- 6 

3- 5 


0.71 (1) 0.13 
0.64 (2) 0.03 


1.21 (4) 8.16 
1.02 (6) 1.62 



On the other hand, comparing the ratio of the fitted "masses" m n /m p for 
the dynamical and the quenched cases (0.92 vs. 0.63) we see that the dynam- 
ical ratio is much higher. In the quenched case such a behavior is observed 
for increased valence-quark masses. Since the dimensionless valence-quark 
mass had the same value 0.1 in both cases we can use the meson mass ratio 
to try to obtain an effective quenched /3i (or lattice spacing) for our dynam- 
ical simulation. Results from the BGR-collaboration's jTT] quenched studies 
at values up to /3i = 8.7 (lattice spacing 0.078 fm) leads us to crudely esti- 
mate the gauge coupling to lie beyond 8.7, corresponding to a lattice spacing 
a < 0.08 fm. 

3.4 Polyakov loop 

For dynamical quarks the Polyakov-loop is not an order parameter and for 
intermediate values of the dynamical quark mass there may be no phase 
transition at all but just a crossover region, analytically connecting the con- 
finement region with the plasma phase. Nevertheless, the Polyakov-loop is 
at least an indicator for the location of that crossover region. In Fig. 0we 
compare the values obtained for the quenched and the dynamical case. We 
find that the center symmetry of the confinement phase is broken for the 
dynamical situation. This confirms our argument that the effective lattice 
spacing has decreased considerably. 

Since we are working on a symmetric lattice of size 8 4 we are not really 
entitled to call T = 1/(8 a) a temperature in a strict sense. However, assum- 
ing a transition temperature range of T c pa 250 MeV (quenched calculations 
give T c pa 270 MeV, dynamical simulations at smaller quark mass give T c pa 
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Figure 7: Scatter plot of the Polyakov loop for the quenched (open symbols) 
and the dynamical (full symbols, m sea = 0.1) situation for lattice size 8 4 and 
A = 7.6. 



170-190 MeV and still comparing it with T we find that a < 0.1 fm 
for our dynamical situation. This is significantly smaller than the quenched 
value 0.19(1) fm. 

This explains the missing zero modes. The dynamical system is too small 
to accommodate instantons and is already in the plasma-like regime. 

4 Conclusions 

We conclude that introducing dynamical fermions for chirally improved ferm- 
ions is possible with reasonable effort. HMC in the simplified version appears 
to work, and dynamical fermions, although still with relatively large mass, 
make a difference in the results as compared to the quenched case. 

For the gauge coupling used we find evidence that the effective lattice 
spacing decreases considerably when switching on dynamical fermions. Al- 
though with our data we cannot derive values for the lattice spacing we have 
several observations indicating the drastic decrease: 

• The Dirac operator spectrum becomes smoother and closer to the GW- 
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circle. Tunneling from non-trivial topological sectors to the trivial one 
is observed, but the system then stays in the sector without zero modes. 

• The mass ratio of pion over rho as derived approximately from the 
correlation functions becomes larger. Since the dimensionless valence- 
quark mass is kept constant this corresponds to a larger physical valence- 
quark mass but smaller lattice spacing, similar to the observations in 
the quenched system. 

• The Polaykov loop shows breaking of the center symmetry. 

All these effects are observed when increasing (3\ in a quenched simulations. 
Curves of constant physics in the m sea ) plane bend towards smaller gauge 
coupling for decreasing m sea . We estimate that the effective lattice spacing 
has changed from 0.19 fm for the quenched simulation at f3\ = 7.6 to a value 
below 0.1 fm for m sea = 0.1. 

Obviously we should work on lattices with larger time extension for better 
analysis of the propagators and at smaller (3\ as a next step. Also desirable 
is an estimate of the growth of the computational effort with decreasing 
sea-quark mass. In view of our results and within these caveats we can be 
optimistic to apply the chiral improved fermion action to a realistic simula- 
tion of QCD. 
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